2+3
a<-5
b<-a+7
a <- 5
a <- 5
> b<-a+7
library("abind", lib.loc="/Library/Frameworks/R.framework/Versions/2.14/Resources/library")
a <- seq(0,50,1)
b<-a+7 + sin(a)
a <- seq(0,50,1)
b <-a+7 + sin(a)
plot(a,b)
plot(a,b, type = 'l')
a <- seq(0,50,1)
b <-a+7 + sin(a^2)
plot(a,b, type = 'l')
a <- seq(0,50,1)
b <-a^2+7 + sin(a)
plot(a,b, type = 'l')
a <- seq(0,50,1)
b <-a+7 + sin(a)
plot(a,b, type = 'l')
boxplot(b)
boxplot(b, col = 'firebrickred2')
boxplot(b, col = 'firebrickred')
boxplot(b, col = 'red')
boxplot(b, col = 'firebrickred')
boxplot(b, col = 'darkseagreen4')
#############
# 1 Objects
#
# R can keep several objects in memory at the same time
# To distinguish them, object have names.
# Objects are assigned wit <-
# Assigning values to objects is done with <-
a <- 5
b <- 6
c <- a * b
a <- 5
b <- 6
c <- a * b
a <- 5
a
b
c
a <- c(1,2,3)
a
cntry <- c("DE", "FR", "NL")
cntry
ls()
rm(cntry)
rm(list=ls())
# scalar
a <- 43
b <- a + 7
a
b
x <- c(a,b,a,b)
y <- x + 10
cntry <- c('Brasil', 'Canada', 'China', 'Singapore')
x
y
cntry
z1 <-cbind(x,y)
z1
z2 <-rbind(x,y)
z2
dat1 <- data.frame(cntry, z1)
dat1
a <- c(1,2,3,4,5)
b <- a + 10
a
a[2]
a[2:4]
m <- cbind(a,b)
m
m[1,1]
m[3,2]
name<- c('mark','luis','peter')
bike<- c('MTB','Single Speed','Race')
name<- c(4,7,8)
name<- c('mark','luis','peter')
bike<- c('MTB','Single Speed','Race')
time<- c(4,7,8)
dat2 <- data.frame(name, bike, time )
dat2
dat2[1,2]
dat2[2,3]
dat2[3,]
x <- c(1,56,23,89,-3,5)
y <- c(24,78,32,27,8,1)
x >20
x[x >20]    # greater as 20
x[x >20 & x !=89] # greater as 20 and unequal 89
x[x>0 | x==-3] # x where x greater 0 oder x=-3
x>0 | x==-3
x==1
y[x==1] # y where x=1
x[x>5 | x < 0]
x[x>5 | x < 2]
##################################################################
##         							                            ##
##  Estimation of presidential positions in LA                  ##
##	 Master File - one file to rule them all...
##  on 2013-08-01				                                ##
##  Copyright (c) 2012 Christian Arnold. All rights reserved.	##
##################################################################
setwd('/Users/Haui/Dropbox/Positions_in_LA/2_models')
# ===========
# = Library =
# ===========
library(foreign)
library(dummies)
library(tm)
library(austin)
# All Bayes Stuff
# Observation: Need to install JAGS in order to run all the bayes packages...
# On a MAC, this means also X11
library(rjags)
library(R2jags)
library(mcmcplots)
library(R2WinBUGS)
# parallel computing
library(snow)
library(dclone)
# Correspondence Analysis
library(ca)
# Plotting
library(ggplot2)
library(gplots)
library(grid)
library(xtable)
# =============
# = Functions =
# =============
# ================================
# = Create Word Frequency Matrix =
# ================================
### Reading in speeches, cleaning the speeches, stemming and then counting
## Function which does all that
#source("code_textprep/preparing_speeches_function.R", chdir = TRUE)
## Implementing the function and writing the word frequency matrices as .csv
#source("code_textprep/preparing_speeches_file.R", chdir = TRUE)
## Load the Frequency Matrices
# Also assigns the tdmat_XXXX and the names.XXX variables
source("code_textprep/data_management_textfiles.R", chdir = TRUE)
#########################################
# Some descriptives on the text corpus ##
#########################################
#source("code_textprep/summarise_corpus.R", chdir = TRUE)
# ==============
# = Estimation =
# ==============
#### Functions
# Wordfish Function from Oli and Jon
#source("estim_positions/functions_wordfish/wordfish.R", chdir = TRUE)
# Function to handle estimation output and convert it to nice data frame
source("estim_positions/functions_wordfish/data_wordfish_out.R", chdir = TRUE)
# Function with filter to use only the most frequent words
# Various thresholds can be implemented: wfm.filter(dat.wfm, keep)
source("code_textprep/wfm_filter.R", chdir = TRUE)
# Filter for frequent words
#source("estim_positions/relative_positions/filter60.R", chdir = TRUE)
#source("estim_positions/relative_positions/filter70.R", chdir = TRUE)
source("estim_positions/relative_positions/filter80.R", chdir = TRUE)
##### Identification via anchors
## Masterfile for retrieving the idealpoints with McKelvey filter from Salamanca Data
## Functions for AMK
# source("estim_positions/Aldrich_McKelvey/amk_functions.R", chdir = TRUE)
## Reads data and calculates idealpoints
# source("estim_positions/Aldrich_McKelvey/amk_master.R", chdir = TRUE)
## This is the general wordfish that takes any number of fixed positions
# source("estim_positions/functions_wordfish/wordfish_fixdoc_ca.R", chdir = TRUE)
## Calculates presidential positions using the AMK filtered positions
# source("estim_positions/relative_positions/idealpoints_and_errors_variance_correction.R", chdir = TRUE)
##### Identification with mean = 0 and sd = 1
##### Estimation via Austin Package
## Estimate with all words
#source("estim_positions/relative_positions/idealpoints_and_errors_austin.R", chdir = TRUE)
## Estimate with 80% most frequent words only
# Estimation
#source("estim_positions/relative_positions/idealpoints_and_errors_austin80.R", chdir = TRUE)
## Estimate with 70% most frequent words only
# Estimation
#source("estim_positions/relative_positions/idealpoints_and_errors_austin70.R", chdir = TRUE)
## Estimate with 60% most frequent words only
# Estimation
#source("estim_positions/relative_positions/idealpoints_and_errors_austin60.R", chdir = TRUE)
# load all estimation results from the RData Sets here
source("estim_positions/relative_positions/load_all_pos_data.R", chdir = TRUE)
### Some descriptives and analyses
# Plots reduced data sets (100%, 80%, 70%, 60%) in one plot
#source("estim_positions/relative_positions/reliability/rel_reducewords.R", chdir = TRUE)
#source("estim_positions/relative_positions/wordtable.R", chdir = TRUE)
################### Bayesian Version ##################
##==============================
## 1 Functions for Estimation
## 1.1 One-core implementation
## Implement until converged
# The JAGS Functions
source("estim_positions/relative_positions/bayes/jags_functions.R", chdir = TRUE)
# The Implementation for 70%
#source("estim_positions/relative_positions/bayes/idealpoints_and_errors_bayes70.R", chdir = TRUE)
# The Implementation for 80%
#source("estim_positions/relative_positions/bayes/idealpoints_and_errors_bayes80.R", chdir = TRUE)
## 1.2 Multicore Implementation
## Chris's office mac runs on an i5, the laptop on an i7 which have 4 and 8 (virtual) cores respectively
## Hence: 'embarassingly parallel' run of 3 chains, implementation on 3 cores
# The JAGS Functions
source("estim_positions/relative_positions/bayes/ep_jags_functions.R", chdir = TRUE)
per.nb.b.80 <- ep.wf.negbin.b(tdmat_per80,nr.chains = 5, nr.iter = 25000)
save(per.nb.b.80, file = '../../../../1_data/dv_prespos/pres_pos_per_negbin_b_80.RData')
per.nb.b.80
setwd('/Users/Haui/Dropbox/Positions_in_LA/2_models/estim_positions/relative_positions/bayes')
save(per.nb.b.80, file = '../../../../1_data/dv_prespos/pres_pos_per_negbin_b_80.RData')
per.nb.b.80.mcmc <- as.mcmc.list(lapply(per.nb.b.80, as.mcmc))
dat <- per.nb.b.80.mcmc
denplot(dat, parms = c("theta[1]"))
per.nb.b.80.mcmc[[2]] <- -1*per.nb.b.80.mcmc[[2]]
denplot(dat, parms = c("theta[1]"))
dat <- per.nb.b.80.mcmc
denplot(dat, parms = c("theta[1]"))
per.nb.b.80.mcmc[[4]] <- -1*per.nb.b.80.mcmc[[4]]
dat <- per.nb.b.80.mcmc
denplot(dat, parms = c("theta[1]"))
per.nb.b.80.mcmc[[4]] <- -1*per.nb.b.80.mcmc[[4]]
dat <- per.nb.b.80.mcmc
per.nb.b.80.mcmc[[3]] <- -1*per.nb.b.80.mcmc[[3]]
dat <- per.nb.b.80.mcmc
denplot(dat, parms = c("theta[1]"))
mcmcplot(dat, parms = c("theta"))
gelman.diag(dat, multivariate = FALSE)
geweke.diag(dat)
#########################################################################
##  Functions for combining the bayes estimates with the existing data ##
##  on 2014-03-6	   			                                       ##
##  Copyright (c) 2014 Christian Arnold. All rights reserved.	       ##
#########################################################################
# Credible Intervals
lowerbound <- function(dat){
quantile(dat, .025)
}
upperbound <- function(dat){
quantile(dat, .975)
}
dat.a.bayes.out <- function(mcmcdat, datout){
# Gets all the indices for the correct thetas on the base of the name
# Collapse multiple chains
mcmcdat <- as.matrix(mcmcdat) 	# This works to collapse multiple chains. Check for a single chain!
allthetas <- grepl('theta',attr(mcmcdat, "dimnames")[[2]])
# this is a function that takes a vector (from apply) and calculates zscores for each element on the basis of the vector
z.vec <- function(data){
(data-mean(data))/sd(data)
}
# This apply takes the data and applies the function z.vec to each row
z.theta <- t(apply(mcmcdat[,allthetas], 1, z.vec))
theta <- mcmcdat[,allthetas]
## Output: Collapse values
# thetas themselves
datout$position.b <- rev(apply(theta, 2, mean))
datout$se.up.b <- rev(apply(theta, 2, upperbound))
datout$se.lb.b <- rev(apply(theta, 2, lowerbound))
datout$move.b <- rev(c(NA,diff(apply(theta, 2, mean))))
datout$move.b[datout$new.pres == TRUE] <- NA			# Correct for the new presidents
datout$abs.move.b <- abs(datout$move.b)
# z scores for thetas
datout$position.b.z <- rev(apply(z.theta, 2, mean))
datout$se.up.b.z <- rev(apply(z.theta, 2, upperbound))
datout$se.lb.b.z <- rev(apply(z.theta, 2, lowerbound))
datout$move.b.z <- rev(c(NA,diff(apply(z.theta, 2, mean))))
datout$move.b.z[datout$new.pres == TRUE] <- NA			# Correct for the new presidents
datout$abs.move.b.z <- abs(datout$move.b.z)
return(datout)
}
dat.dimcorrection.b <- function(dat){
# correct values
dat$position.b 	<- -1*dat$position.b
dat$se.up.b 		<- -1*dat$se.up.b
dat$se.lb.b 		<- -1*dat$se.lb.b
dat$move.b 		<- -1*dat$move.b
# correct z scores
dat$position.b.z 	<- -1*dat$position.b.z
dat$se.up.b.z 		<- -1*dat$se.up.b.z
dat$se.lb.b.z 		<- -1*dat$se.lb.b.z
dat$move.b.z 		<- -1*dat$move.b.z
return(dat)
}
dat.pres.pos.per80 <- dat.a.bayes.out(per.nb.b.80.mcmc, dat.pres.pos.per80)
dat.pres.pos.per80
save(dat.pres.pos.per80, file = '../../../../1_data/dv_prespos/pres_pos_per80.RData')
dat.pres.pos.all80 <- rbind(dat.pres.pos.arg80
,dat.pres.pos.bra80
,dat.pres.pos.chi80
,dat.pres.pos.col80
,dat.pres.pos.cri80
,dat.pres.pos.ecu80
,dat.pres.pos.gtm80
,dat.pres.pos.mex80
,dat.pres.pos.per80
,dat.pres.pos.pry80
,dat.pres.pos.slv80
,dat.pres.pos.ury80
,dat.pres.pos.ven80)
save(dat.pres.pos.all80, file = '../../../../1_data/dv_prespos/pres_pos_all80.RData')
write.dta(dat.pres.pos.all80, file = "../../../../1_data/ADW_MasterData.dta")
